import holoviews as hv
hv.extension('bokeh')
import numpy as np
import scipy.signal
from IPython.display import YouTubeVideo,  Audio

1. Introducción al procesamiento digital de señales

1.1. Sistemas y señales

Definición de sistema: Es un proceso que produce una señal de salida en función de una señal de entrada

../../_images/signal-intro.png

Definición de señal: Es una función que describe o representa el comportamiento de un fenómeno físico o sistema. Este comportamiento puede variar ya sea en el tiempo, espacio u otra variable

1.1.1. Ejemplos de señales

Revisemos los siguientes ejemplos ¿Cuáles son las variables independientes y dependientes en cada caso?

Las variaciones del IPSA entre 2006 y 2015

../../_images/signal-ipsa.png

Los pulsos de reloj de un microcontrolador

../../_images/signal-clock.png

El brillo relativo de una estrella medido desde la Tierra. La estrella es orbitada por un exoplaneta

../../_images/signal-transit.gif

Una conversación humana grabada digitalmente

../../_images/signal-audio.gif

La actividad cerebral de una persona medida usando ElectroEncefalografía

../../_images/signal-eeg.png

Datos meteorológicos de un huracán

../../_images/signal-weather.gif

Una película

../../_images/signal-bunny.gif

1.1.2. ¿Qué es el procesamiento de señales?

Es la disciplina que se dedica al diseño de sistemas para

  • representar …

  • filtrar …

  • codificar …

  • transmitir …

  • estimar …

  • detectar …

  • inferir …

  • descubrir …

  • reconocer …

  • síntesizar …

  • reproducir …

  • etc

… señales

Fuente: Jose S.F. Moura, “What is signal processing”, IEEE signal processing magazine

YouTubeVideo('R90ciUoxcJU')

1.2. Definición matemática de una señal

Una señal se representa matematicamente como una función

\[ y = f(x) \]

es decir un mapeo entre el espacio de entrada y el espacio de salida

\[ f: \mathcal{X} \rightarrow \mathcal{Y}, \quad x \in \mathcal{X} \wedge y \in \mathcal{Y}, \]

donde:

  • \(x\) se llama variable independiente, entrada o argumento y su espacio se llama dominio

  • \(y\) se llama variable dependiente, salida o retorno y su espacio se llama codominio

1.2.1. Funciones matemáticas típicas

Algunas funciones típicas que veremos a lo largo de este curso son

  • Impulso unitario:

\[\begin{split} \delta[n] = \begin{cases} 1 & n=0\\ 0 & n \neq 0 \end{cases} \end{split}\]
  • Escalón unitario:

\[\begin{split} u[n] = \sum_{k=0}^\infty \delta[n-k] = \begin{cases} 1 & n\geq 0\\ 0 & n < 0 \end{cases} \end{split}\]

Propiedad: \(\delta[n] = u[n] - u[n-1]\)

  • Exponencial real:

\[ y[n] = C e^{\alpha n} \]
  • Exponencial compleja:

\[ y[n] = e^{jn\omega} = \Re [e^{jn\omega}] + j \Im[e^{jn\omega}] = \cos(n\omega) + j \sin (n\omega) \]
  • Gaussiana

\[ y[n] = e^{- \alpha (n-n_0)^2}, \alpha> 0 \]
x = np.linspace(-4, 4, num=1000)
p1 = hv.Overlay([hv.Curve((x, x), label='y = x'), hv.Curve((x, np.abs(x)), label='y = |x|')])
p2 = hv.Overlay([hv.Curve((x[x>0], np.sqrt(x[x>0])), label='y = sqrt(x)'), 
                 hv.Curve((x[x>0], np.log(x[x>0])), label='y = log(x)')])
p3 = hv.Overlay([hv.Curve((x, np.exp(x/4)), label='y = exp(x/4)'), 
                 hv.Curve((x, np.exp(-x/4)), label='y = exp(-x/4)')])
p4 = hv.Overlay([hv.Curve((x, np.cos(x)), label='y = cos(x)'), 
                 hv.Curve((x, np.sin(x)), label='y = sin(x)')])
p5 = hv.Overlay([hv.Curve((x, np.heaviside(x, 0)), label='y = step(x)'), 
                 hv.Curve((x, 1.0/(1.0+np.exp(-x))), label='y = logit(x)')])
p6 = hv.Overlay([hv.Curve((x, scipy.signal.unit_impulse(len(x), 'mid')), label='y = delta(x)'), 
                 hv.Curve((x, np.exp(-0.5*x**2)), label='y=exp(-x^2)')])
plot = hv.Layout([p1 + p2 + p3] + [p4 + p5 + p6])
plot.opts(hv.opts.Curve(width=280, height=280, xlim=(-4,4), ylim=(-4,4), show_grid=True, tools=['hover']), 
          hv.opts.Overlay(legend_position='top')).cols(3)

1.2.2. Algunas propiedades de las funciones

  • Función par o simétrica:

\[ f(x) = f(-x) \]
  • Función inpar o antisimétrica:

\[ f(x) = -f(-x) \]
  • Función periódica: Sea \(T\) el período de la función

\[ f(x) = f(x+T) \]
  • Función lineal: Sean \(a\), \(b\) coeficientes escalares

\[ f(ax + by) = a f(x) + b f(y), \]

El periódo de una función es equivalente al inverso de su frecuencia fundamental

La siguiente animación muestra como varía una sinusoide en función de su frecuencia

time = np.linspace(-6, 6, num=1000)

hmap = hv.HoloMap(kdims=['Frecuencia'])
for freq in np.arange(0, 1, step=0.01):
    hmap[freq] = hv.Curve((time, np.sin(2*np.pi*freq*time)))

hmap.opts(hv.opts.Curve(width=500, height=300))
hv.output(hmap, holomap='gif', fps=10)

1.3. Señales analógicas y digitales

La siguiente figura muestra cuatro combinaciones en que se pueden presentar las variables independiente y dependiente de una señal

  • Variable independiente continua (columna izquierda)

  • Variable independiente discreta (columna derecha)

  • Variable dependiente continua o análoga (fila superior)

  • Variable dependiente discreta o cuantizada (fila inferior)

../../_images/signal-classification1.png

Las señales “naturales” son en general analógicas de tiempo continuo. Sin embargo los computadores trabajan con señales digitales de tiempo discreto

Podemos convertir una señal analógica a una representación digital con los siguientes pasos:

  • Se discretiza en el tiempo muestreando según el tiempo de muestreo \(T_s\) del sistema, por ejemplo:

\[ x_A(t=kT_S) = x[k], \quad k \in \mathbb{Z} \]

Esto convierte la señal analógica \(x_A\) en una secuencia discreta de valores \(x\) indexada por un entero \(k\), es decir un arreglo

  • Se cuantiza la amplitud de la señal, por ejemplo:

\[\begin{split} v_D = \begin{cases} 0 & v_A \in [0.0, 0.8] V \\ 1 & v_A \in [2.0, 5.0] V\end{cases} \end{split}\]

Esto convierte la variable dependiente continua (en este caso voltaje) a una representación discreta (en este caso binaria)

Ts = 0.249
x = np.linspace(0.0, 2.0, num=1000)
y = np.sin(2.0*np.pi*1.0*x)

p1 = hv.Curve((x, y), label='Señal analógica').opts(width=250)
p2 = hv.Curve((x, Ts*np.round(y//Ts)), label='Señal cuantizada').opts(width=250)
p3 = hv.Points((x[::20], Ts*np.round(y[::20]//Ts)), label='Señal digital').opts(width=250)
(p1 + p2 + p3)

1.4. Propiedades fundamentales de las señales

En general trabajaremos con señales que, dentro de un rango de interés, son acotadas en energía, potencia y/o en ancho de banda. Es decir son señales finitas que no divergen ni se vuelven singulares

La energía de una señal mide su “tamaño” o el “espacio que ocupa”. Para una señal analógica y discreta, respectivamente

Matemáticamente, se define como

\[ E_s = \int_{-\infty}^\infty |s(t)|^2 \,dt \qquad E_s = \sum_{n=-\infty}^\infty |s[n]|^2 \]

Una señal acotada en energía debe cumplir \(E_s < \infty\)

La potencia promedio de una señal se define como su energía por unidad de tiempo

Matemáticamente, se define como

\[ P_s = \lim_{T\to \infty} \frac{1}{2T} \int_{- T}^{T} |s(t)|^2 \,dt \qquad P_s = \lim_{N\to \infty} \frac{1}{2N+1} \sum_{n = - N}^{N} |s[n]|^2 \]

Una señal acotada en potencia debe cumplir \(P_s < \infty\) y una señal de energía finita tiene potencia cero

Por otro lado una señal de potencia finita tiene duración infinita

El ancho de banda de una señal mide su tasa de cambio o velocidad

Una señal acotada en ancho de banda debe tener transiciones suaves

Ejemplo formativo

¿Cual es la energía de esta señal?

\[\begin{split} s(t) = \begin{cases} 0 & t < 0 \\ 2e^{-t/2} & t \geq 0\end{cases} \end{split}\]
dt = 0.001
x = np.arange(-5, 20, step=dt)
y = np.zeros(shape=(x.shape))
y[x>=0] = 2*np.exp(-x[x>=0]/2)
# Integral numérica
print((y**2).sum()*dt)

hv.Curve((x, y)).opts(width=500)
4.0020003250765726

¿Puedes comprobar el resultado resolviendo la integral?

1.5. Procesamiento estadístico de señales

Las funciones matemáticas que vimos anteriormente son ejemplos de señales deterministas. Una señal es determinista si puede describirse completamente por una ecuación matemática

Por ejemplo

\[ y = \cos(x)\sin(2x) + 0.1x \]
x = np.linspace(-10, 10, num=1000)
y = np.cos(x)*np.sin(2*x) + 0.1*x
hv.Curve((x, y)).opts(width=500)

Sin embargo, existen algunos fenómenos que son difíciles o imposibles de describir de forma determinísta, como son por ejemplo las señales estocásticas

Una señal estocástica tiene un cierto grado de aleatoridad, para describirla debemos usar teoría de probabilidades

Por ejemplo

\[ z \sim \mathcal{N}(\mu, \Sigma) \]

quiere decir que \(z\) se distribuye normal con media \(\mu\) y covarianza \(\Sigma\)

Podemos estudiar la distribución de \(z\) usando un histograma

z = np.random.multivariate_normal(np.zeros(100), np.eye(100))
edges, bins = np.histogram(z, bins=20)
hv.Histogram((edges, bins)).opts(width=500)

1.5.1. Breve repaso de estadística

Variable aleatoria (VA): Es una variable cuyos valores posibles son resultados de un fenomeno aleatorio

  • Se describe en términos de su distribución

    • función de densidad de probabilidad (fdp) para variables continuas

    • función de masa de probabilidad (fmp) para variables discretas

  • Los valores observados a partir de la VA se llaman realizaciones

  • Por lo general asumimos que las realizaciones son iid: independientes e identicamente distribuidas

  • Usualmente se denota como \(X\) (mayúscula) y sus realizaciones como \(\{x_1, x_2, \ldots, x_N\}\) (minúscula)

Una distribución se describe a través de sus momentos estadísticos. Para una variable aleatoria \(X\) su momento de orden \(k\) es

\[ \mu_k = \mathbb{E} \left[X^k\right] \]

y su momento central de orden \(k\)

\[ \hat \mu_k =\mathbb{E} \left[(X - \mathbb{E}[X])^k\right] \]

donde la esperanza se define como

\[ \mathbb{E}[X] = \sum_n x_n p(x_n) \]

y la media muestreal es

\[ \langle X \rangle = \frac{1}{N} \sum_{n=1}^N x_n \]

Los momentos estadísticos más usados son la

  • Media (\(\mu_1\)): Define la localización de la distribución

  • Varianza (\(\hat \mu_2\)): Define la dispersión o ancho de la distribución

  • Simetría (\(\frac{\hat \mu_3 }{\sqrt{\hat \mu_2^3}}\)): Define el peso relativo de las colas de la distribución

  • Kurtosis (\(\frac{\hat \mu_4 }{ \hat \mu_2^2}\)): Define el peso relativo entre las colas y el centro (moda) de la distribución

La siguiente gráfica muestra como cambia una distribución ante estos momentos

../../_images/moments.png

1.5.2. Proceso aleatorio/estocástico

Un proceso aleatorio (PA) es una colección de variables aleatorias indexadas

\[ X^t = {X_1, X_2, X_3, \ldots, X_T} \]

Llamamos a una realización del PA: serie de tiempo

\[ X^t \sim \{x_t\}_{t=1,\ldots,T} = {x_1, x_2, x_3, \ldots, x_T} \]

Un PA se dice estacionario si sus momentos no varían con el tiempo, es decir

\[ \mathbb{E}[X_1] = \mathbb{E}[X_2] = \ldots = \mathbb{E}[X_T] \]

Un PA se dice ergódico si los momentos se pueden deducir a partir de una realización (suficientemente larga) del mismo

  • Esto es útil pues muchas veces observamos sólo una realización

  • Los promedios muestreales del ensamble equivalen a promedios muestreales temporales

\[ \langle x_t \rangle = \mathbb{E}[X_k] ~~\forall k \]

La siguiente figura muestra la diferencia entre un promedio temporal y un promedio de las realizaciones

../../_images/stationary-and-ergodic.png

Ejemplo: Realizaciones de un PA con distribución Gaussiana

Mientras más muestras observamos, mejor se aprecia la distribución subyacente

np.random.seed(12345)
data = np.random.randn(10000)

x = np.linspace(-3, 3, num=100); 
signal = hv.HoloMap(kdims=['Muestras'])
histogram = hv.HoloMap(kdims=['Muestras'])
for N in [10, 50, 100, 500, 1000, 5000, 10000]:
    signal[N] = hv.Curve((data[:N]), kdims='tiempo', vdims='datos').opts(xlim=(1, N+1), logx=True)
    edges, bins = np.histogram(data[:N], bins=20, range=(-3, 3), density=True)
    histogram[N] = hv.Overlay([hv.Histogram((edges, bins), kdims='datos', vdims='Frecuencia').opts(alpha=0.5),
                               hv.Curve((x, np.exp(-x**2)/np.sqrt(2.0*np.pi)))])
    
hv.output((signal + histogram), holomap='gif', fps=1)

Ejemplo: Realizaciones de un PA con distribución Normal/Gaussiana con media cero

Cada columna tiene una covarianza distinta. Mientras más grande es el valor de “tau” más correlacionadass están las muestras

np.random.seed(12345)
dt = np.repeat(np.reshape(np.arange(100), (1, -1)), 100, axis=0)

def draw_random_plot(tau):
    p = []
    for k in range(0, 25, 5): 
        Sigma = np.exp(-0.5*np.square(dt - dt.T)/tau**2)
        data = k + np.random.multivariate_normal(np.zeros(100), Sigma)
        p.append(hv.Curve((data)))        
    return hv.Overlay(p, label=f'tau = {tau}').opts(width=250, height=400, yaxis=None, xaxis=None)

draw_random_plot(0.5) + draw_random_plot(1) + draw_random_plot(2)

1.5.3. Ruido

En la práctica no solemos observar señales puramente deterministas. Una de las causas de esto es la presencia del ruido

El ruido es una corrupción indeseable y usualmente estocástica que modifica la señal de interés

Corrupción aditiva

Un supuesto muy utilizado es considerar que el ruido es una corrupción aditiva

\[ y[k] = x[k] + n[k], \]

donde \(y\) es la señal observada, \(x\) es la señal de interés y \(n\) es una señal de ruido

np.random.seed(12345)
t = np.linspace(-10, 10, num=500)
x = np.cos(t)*np.sin(2*t) + 0.1*t
s = 0.5
n = s*np.random.randn(len(t))

print(f"SNR: {10*np.log10(np.sum(x**2)/np.sum(n**2))}")

p1 = hv.Curve((t, x))
p2 = hv.Curve((t, n))
p3 = hv.Curve((t, x+n))
(p1+p2+p3).opts(hv.opts.Curve(width=250, height=250))
SNR: 4.063138926337416

Podemos cuantificar el “nivel de ruido” o la claridad de la señal observado en términos de la razón señal a ruido (Signal to Noise Ratio, SNR)

  • La SNR se define como la razón entre la energía de la señal y la energía del ruido

  • La SNR se mide en decibeles [dB]

\[ \text{SNR} = 10 \log_{10} \frac{E_x}{E_n} \]

1.5.4. Tipos de ruido

Podemos clasificar el ruido según

  • la distribución que sigue, por ejemplo ruido Gaussiano o ruido Uniforme

  • la dependencia temporal/espectral entre sus realizaciones, por ejemplo ruido blanco o ruido rojo

¿Cómo se ve y suena un ruido blanco?

Fs = 44100
y = np.random.randn(Fs*5)
y = y/np.max(np.absolute(y))

hv.Curve((y), 'tiempo', 'ruido').opts(width=500)

¿y un ruido rojo?

dy = np.random.uniform(-1, 1, Fs*5)
y = np.cumsum(dy)
y = (y - np.mean(y))
y = y/np.amax(np.absolute(y))
hv.Curve((y), 'tiempo', 'ruido').opts(width=500)
Audio(y, rate=Fs)

Volveremos a estudiar sobre el “color” del ruido cuando estudiemos la Transformada de Fourier

1.5.5. Ejemplo formativo: Reducción de ruido usando promedios

Sea una señal periódica con ruido aditivo, estacionario y de media cero

En este caso podemos reducir el ruido haciendo promedios de la señal de interés

\[ \mathbb{E}[y] = \mathbb{E}[x + n] = \mathbb{E}[x] + \mathbb{E}[n] = \mathbb{E}[x] \]

Notar que esto sólo funciona si alineamos adecuadamente la señal según su período. Además el ruido sólo se hace cero cuando el número de muestras promediadas tiende a infinito

Fs = 200 # Frecuencia de muestro
time = np.arange(0, 50, step=1/Fs) 
T, s = 1., 2.
# La señal de ejemplo es una onda cuadrada con ruido blanco gaussiano
data = scipy.signal.square(2.0*np.pi*time/T) + s*np.random.randn(len(time));     
hv.Curve((time, data), 't', 'x').opts(width=500)

Si cortamos la señal en trozos de tamaño \(T\) y tomamos el promedio obtenemos

hv.Curve((data.reshape(-1, int(Fs*T)).mean(axis=0)), 't', 'y')

Se aprecia claramente la forma cuadrada de la señal.

A lo largo de este curso veremos métodos y algoritmos más robustos para filtrar y procesar señales contaminadas con ruido